Insights into the phylogeny of the ciliate of class Colpodea based on multigene data

Abstract In the class Colpodea, there are many unresolved evolutionary relationships among taxa. Here, we report 30 new sequences including SSU‐rRNA, ITS1‐5.8S‐ ITS2 rRNA, and the mitochondrial small subunit ribosomal RNA (mtSSU‐rRNA) genes of five colpodeans, and conduct phylogenetic analyses based on each individual gene and a two‐gene concatenated dataset. For the first time, multi‐genes were used to analyze phylogenetic relationships in the class Colpodea. The main findings are: (1) SSU‐rRNA, ITS1‐5.8S‐ ITS2 rRNA, and mtSSU‐rRNA gene sequences of C. reniformis and C. grandis are provided for the first time, and these two species group into the clade including C. inflata, C. lucida, C. cucullus, and C. henneguyi; (2) clustering pattern and morphological similarity indicate that Bresslauides discoideus has a close relation with Colpodidae spp.; (3) Emarginatophrya genus diagnosis is improved to be ‘Hausmanniellidae with sharply shortened and isometric leftmost 1‐4 ciliary rows’ and Colpoda elliotti is transferred to Emarginatophrya; (4) the genus Colpoda is still non‐monophyletic with the addition of 10 populations from five Colpoda species sequences, but there are only two Colpoda groups left based on the present work: Group I comprises C. inflata, C. lucida, C. cucullus, C. henneguyi, C. reniformis, and C. grandis, Group II comprises C. maupasi and C. ecaudata, and the presence of diagonal grooves and the way the vestibular opens might be the two key features that differentiates Colpoda species groups; (5) a close molecular relationship, and highly similar merotelokinetal mode, somatic ciliary pattern, and basic organization of the oral apparatus with P. steinii suggests Bromeliothrix metopoides should be temporarily assigned to Colpodidae.


| INTRODUC TI ON
Ciliates are a highly differentiated group of microbial eukaryotes that inhabit virtually all environments on the Earth's surface where there is sufficient water for their survival. They play a fundamental role in microbial food web function by mediating the transfer of organic matter and energy between different trophic levels in a wide range of ecosystems (Agatha et al., 2021;Bai et al., 2020;Chi et al., 2021;Dias et al., 2021;Liu et al., 2022;Song et al., 2021;Wang et al., 2021;Weisse & Montagnes, 2021;Wu et al., 2021). Due to their particular nature (e.g., high reproductive capacity, diverse morphology, nuclear dimorphism, chromosomal fragmentation, and sensitivity to environmental changes), ciliates are important model organisms in many disciplines, such as epigenetics, symbiotic relationships, molecular biology, organismal development, evolution, biogeography, and ecology (Berger, 2011;Cheng et al., 2019;Corliss, 1979;Gao et al., 2014;Lynn, 2008;Song et al., 2009;Zhao et al., 2020).
Characteristically, colpodeans produce resting cysts that allow them to survive desiccation and other adverse environmental conditions and thus often constitute a significant hidden soil biodiversity component (Foissner, 1987;Quintela-Alonso et al., 2011). Although their oral apparatus is highly diverse, colpodeans share a similar somatic cortex structure which has led to misclassification among species (Foissner, 1993;Lynn, 2008;Quintela-Alonso et al., 2011;Vďačný & Foissner, 2019). Small and Lynn (1981) recognized colpodeans as a monophyletic taxon and established the class Colpodea on the basis of a set of unique morphological characteristics. The application of molecular techniques supported Colpodea monophyly and its position within the subphylum Intramacronucleata Lynn, 2003), but also resulted in a major re-evaluation of interclass classification and evolution Vďačný & Foissner, 2019).
Recently, investigations based on multigene trees have been increasingly used for phylogenetic studies and generally produce more robust results than those based on single gene markers (Gao et al., 2013(Gao et al., , 2016Sun et al., 2020;Zhang et al., 2020). Multigene analyses usually include nSSU-rRNA, ITS1-5.8S-ITS2 and nLSU-rRNA genes, which are in the same chromosome . However, Colpodea phylogenetic analyses so far have focused on a single gene, (i.e., the nSSU-rRNA gene or the mtSSU-rRNA gene), and few studies have used a multigene approach (e.g., Dunthorn et al., 2011). Furthermore, most colpodean sequences in the NCBI database are of SSU-rRNA, while other gene sequences are very limited. This is one of the main dilemmas of colpodean phylogenetics.
With the conservative evolution of SSU rDNA alongside various issues such as asynchronous evolution with morphology, delineation of Colpodea species remains problematic. And, the usefulness of combination of nuclear and mitochondrial gene is rarely tested in Colpodea. The mtSSU-rRNA gene is more variable and has a higher evolutionary rate compared to the nSSU-rRNA gene (Boore & Brown, 1998;Moore, 1995;, and it is reasonable to assume that it can effectively uncover deep nodes within Colpodea Rand, 2003;Zhang et al., 2019). In this study, 30 new sequences of 10 populations from five colopdeans (Colpoda cf. inflata, C. reniformis, C. inflata pop1-2, C.  with rice or wheat grains added to enrich bacteria growth as a food source. Clonal cultures were established for each of the 10 populations from five colpodean species isolated. Silver carbonate staining (Foissner, 1992) was used to reveal the infraciliature.

| DNA extraction, amplification, and sequencing
For each species, 1-10 cells from clonal culture were isolated under a stereomicroscope using a micropipette, washed with distilled water at least thrice to remove potential contaminants, and then incubated in non-nutrient distilled water for 24 h. Cells were then transferred to an Eppendorf tube in a volume of no more than 5 μl distilled water. Total genomic DNA was extracted using the DNeasy & Tissue Kit (Hilden, QIAGEN) following the manufacturer's instructions.
SSU-rRNA, ITS1-5.8S-ITS2 rRNA, and mtSSU-rRNA genes were amplified by the polymerase chain reaction (primers see Appendix Table A1). High-fidelity Taq polymerase (Takara Ex Taq; Takara Biomedicals) was used to reduce amplification errors. PCR condition for the SSU-rRNA gene amplification was denaturation for 5 min at 94°C, followed by 5 cycles of denaturation for 30 s at 94°C, annealing for 1 min 45 s at 56°C, extension for 2 min at 72°C and the other 25 cycles of denaturation for 45 s at 94°C, annealing for 1 min 45 s at 60°C, extension for 2 min at 72°C and a final extension at 72°C for 8 min. The ITS1-5.8S-ITS2 rRNA gene was amplified as follows: 5 min initial denaturation 94°C, followed by 35 cycles of denaturation for 30 s at 94°C, annealing for 45 s at 58°C, extension for 1 min at 72°C and a final extension at 72°C for 10 min. The mtSSU-rRNA gene was amplified as follows: 5 min initial denaturation 94°C, 5 cycles of 45 s at 94°C, 1 min 45 s at 58°C and 2 min at 72°C, and the other 25 cycles of denaturation for 45 s at 94°C, annealing for 1 min 45 s at 60°C, and extension for 2 min at 72°C with a final extension of 10 min at 72°C. Cloning and sequencing were both routine operations (Zhang et al., 2019). A total of 30 new SSU-rRNA, ITS1-5.8S-ITS2 rRNA, and mtSSU-rRNA gene sequences from five Colpodea species were determined (GenBank accession numbers see Appendix Table A2).

| Data sets and alignments
Other colpodeans and outgroup sequences were obtained from GenBank, and their accession numbers are shown in Figures 1-4.
Sequences downloaded with questionable morphological information, unknown source, or of too short length were not included. SSU-rRNA gene was aligned using Clustal W implemented in BioEdit 7.0.1 (Hall, 1999). The ITS1-5.8S-ITS2 rRNA gene and mtSSU-rRNA gene sequences were aligned using the default parameters implemented in Guidance server (http://guida nce.tau.ac.il/) (Penn et al., 2010).
To remove ambiguously aligned positions in the multiple sequence alignments all the sequence datasets were edited by eye in BioEdit 7.0.1 (Hall, 1999). Sequences used for phylogenetic analyses were compiled into four data sets: (i) 1833 characters of SSU-rRNA (83 taxa); (ii) 722 characters of ITS1-5.8S-rRNA -ITS2 (19 taxa); (iii) 1138 characters of mtSSU-rRNA (31 taxa); and (iv) 2895 characters of concatenated sequence data of the above two genes (29 taxa in total).
The two alignments were concatenated to be contiguous in SeaView V4 for gene analysis. All new sequences were deposited in the NCBI database (https://www.ncbi.nlm.nih.gov/), accession numbers, lengths, and G&C contents are shown in Appendix Table A2. For all trees, the outgroup was composed of three species of tetrahymenid taxa (Tetrahymena pyriformis, T. tropicalis and T. americanis).

| Phylogenetic analyses based on SSU-rRNA gene sequence data
Each of the Colpodea four orders are monophyletic and the genus Colpoda is non-monophyletic. The orders Colpodida and Cyrtolophosidida cluster together to form a clade that is a sister to Bursariomorphida. The order Platyophryida is the deepest branching lineage within the Colpodea.

F I G U R E 1
The maximum-likelihood (ML) tree based on the concatenated genes (SSU-rRNA and mtSSU-rRNA genes) of major members of the class Colpodea. Newly added sequences in this study are bolded in red. Node support is shown as: ML bootstraps/BI posterior probability. "-" indicate mismatch in topology between Bayesian and ML tree. Fully supported (100%/1.00) branches are marked with solid circles. The scale bar corresponds to five substitutions per 100 nucleotide sites.

F I G U R E 2
The maximum-likelihood (ML) tree based on the SSU-rRNA gene of major members of the class Colpodea. Newly added sequences in this study are bolded in red type. Node support is shown as: ML bootstraps/BI posterior probability. "-" indicate mismatch in topology between Bayesian and ML tree. Fully supported (100%/1.00) branches are marked with solid circles. Two long branches has been shortened, as shown by "//", and the other branches are drawn to scale. The scale bar corresponds to five substitutions per 100 nucleotide sites. (a-j) The five newly sequenced species after silver carbonate staining. Paracolpoda steinii sequences including four newly sequenced species (P. steinii pop1-4) cluster together with Bromeliothrix metopoides.

| Phylogenetic analyses based on ITS1-5.8S-ITS2 region sequence data and the secondary structures of ITS2
In the ITS1-5.8S-ITS2 rRNA gene tree, bootstrap values for the maximum likelihood (ML) and posterior probabilities for the Bayesian inference (BI) were mapped onto the best BI tree (Figure 3). Due to sequence shortages, only the order Colpodida (family Colpodidae, genus Colpoda and Paracolpoda) and outgroup were included in this tree. Colpodida, Colpoda, and Paracolpoda form a monophyletic clade. Colpoda reniformis is still separated from the other newly sequenced Colpoda species (C. inflata pop1-2 and Colpoda cf. inflata), and forms a sister clade with the newly sequenced C. grandis (96% ML, 1.00 BI). Two populations of C. inflata (KM222071 and MZ562700), C. inflata pop1, and Colpoda cf. inflata cluster together, and the clade then forms a sister clade with C. inflata. Colpoda inflata pop2 cluster with the C. inflata + C. inflata pop1 + Colpoda cf.
inflata clade. Seven Paracolpoda species cluster in one clade, and the clade formed by the five newly sequenced Paracolpoda steinii cluster as a sister group to P. steinii (AB684403) with moderate support (77% ML, 1.00 BI). Paracolpoda steinii (AB684403) groups with the F I G U R E 3 The Bayesian inference (BI) tree based on the ITS1-5.8S-ITS2 gene of major members of the class Colpodea. Newly added sequences in this study are bolded in red type. Node support is shown as: ML bootstraps/BI posterior probability. "-" indicate mismatch in topology between Bayesian and ML tree. Fully supported (100%/1.00) branches are marked with solid circles. The scale bar corresponds to 0.2 expected substitutions per site.

F I G U R E 4
The Bayesian inference (BI) tree based on the mtSSU-rRNA gene of major members of the class Colpodea. Newly added sequences in this study are bolded in red type. Node support is shown as: ML bootstraps/BI posterior probability. "-" indicate mismatch in topology between Bayesian and ML tree. Fully supported (100%/1.00) branches are marked with solid circles. The scale bar corresponds to 0.1 expected substitutions per site.

F I G U R E 5
The putative secondary structures of ITS2 in the present study.

| Phylogenetic analyses based on mtSSU-rRNA gene sequence data
The BI tree topology is similar to the ML tree, so only the BI tree with branch-support values is shown in Figure 4

| Phylogeny of genus Colpoda
In previously reported SSU rRNA gene trees, typical Colpoda species are distributed over the whole Colpodea tree, and mainly form four groups: C. inflata + C. cucullus + C. lucida, C. henneguyi, C. maupasi + C. ecaudata, and C. elliotti clades Foissner & Stoeck, 2009;Lynn et al., 1999;Vďačný & Foissner, 2019). This motivated Foissner et al. (2011) and Dunthorn et al. (2012) to suggest rapid Colpoda radiation that produced several new genera and families. We added 30 new multiple gene sequences of 10 populations from five colopdeans. These additions change the Colpoda species into three main groups: Group I is composed of C. inflata, two Chinese C. inflata populations, Colpoda cf. inflata, C. lucida, C. cucullus, C. henneguyi, C. reniformis, and C. grandis; Group II comprises C. maupasi and C. ecaudata; and Group III only contains C. elliotti in the SSU rRNA gene trees.
Colpoda reniformis and C. grandis clustered with C. henneguyi, and the resulting clade then clustered with Bresslauides discoideus. In this case, B. discoideus groups into the Colpoda clade including C. inflata, C. lucida, C. cucullus, C. henneguyi, C. reniformis, and C. grandis in the SSU rRNA, mtSSU-rRNA gene, and concatenated trees. Bresslauides species differ mainly by the unique right semicircular oral polykinetid that is longer than the left, as opposed to being of equal length in the Colpodidae (Foissner, 1987(Foissner, , 1993, but similar with the features of Colpodidae sp. in vivo, which indicates that Bresslauides species have a close relation with Colpodidae sp. Broader sampling and molecular sequencing of Bresslauides species would be necessary to clarify this relationship, although we also believe that the character of the right semicircular curved oral polykinetid may have evolved more than once (Bourland et al., 2011;Dunthorn et al., 2008).

Colpoda elliotti KJ873047 clusters in the clade formed by
Kreyellidae-like species, Emarginatophrya aspera, and it unites with E. aspera in trees comprised colpodids only . However, establishing this new colpodid genus has not helped to completely erase the Colpoda paraphyly problem. Colpoda aspera and C. elliotti are highly similar morphologically and according to Foissner (1993), 'the left oral polykinetid of C. elliotti is slightly variform'. Then, we found in fig. 53b in Foissner's book that the left oral polykinetid was slightly emarginated (as a result of its leftmost kinety always consisting of only two basal bodies) in some representative individuals. The only two Emarginatophrya species, E. aspera and E. terricola, also have three to four sharply shortened leftmost ciliary rows (almost of equal length), which were described as 'emarginated left oral polykinetid'. Therefore, we decided to change the diagnosis of Emarginatophrya to 'Hausmanniellidae with sharply shortened and isometric leftmost 1-4 ciliary rows'. So, we transferred
As a result, there are now two remaining Colpoda groups: Group I comprises C. inflata, C. lucida, C. cucullus, C. henneguyi, C. reniformis, and C. grandis, and Group II comprises C. maupasi and C. ecaudata. Colpoda maupasi and C. ecaudata group with Exocolpoda and Ropoma species, which have distinctly different characteristics from the former species, including different life cycles, boomerang-shaped left oral ciliary fields, and a very thick resting cyst wall. Colpoda maupasi and C. ecaudata are two common Colpoda species, but they group far away from other congeners. A probable explanation for this grouping pattern is that both C. maupasi and C. ecaudata lack the diagonal groove and their vestibular opening is marked on the left body margin by shallow indentation, but other Colpoda species having molecular information possess diagonal grooves (except C. inflata, but its left body margin with right-angular notch, in the center where the vestibulum opens; Foissner, 1993). Morphologically, Tillina magna also possess diagonal grooves (Foissner, 1993), which might explain why Tillina magna and the Group I of the genus Colpoda cluster together in the SSU-rRNA and mtSSU-rRNA gene trees.
Collectively, the presence of diagonal grooves and the way the vestibular opens might be the two key features that differentiate Colpoda species groups. Due to sequence shortages, only one species of Tillinidae used in the concatenated tree, which leads to the clustering of Tillina magna with Emarginatophrya aspera.
Therefore, we suggest that sampling should be expanded to further address this issue.

| Paracolpoda and Bromeliothrix phylogeny
Totally twelve populations (including five from the present work) of Paracolpoda steinii firmly cluster with Bromeliothrix metopoides in the SSU rRNA gene trees. Bromeliothrix metopoides has the same merotelokinetal mode with Paracolpoda steinii, and its somatic ciliary pattern is similar to that of Colpoda species (Foissner, 1993(Foissner, , 2010Weisse et al., 2013). Bromeliothrix has a ciliary and silverline pattern typical for members of the family Colpodidae, and its basic organization of the oral apparatus is also the same as that known from genera of Colpodidae families (Foissner, 2010;Foissner et al., 2002). Foissner pointed out that Bromeliothrix neither belong to the Exocolpodidae nor the Hausmanniellidae due to its complex morphological characteristics, and several phylogenetic analyses based on SSU rRNA gene sequence also indicated that B. discoideus was more closely related to a clade formed by large Colpoda species than with Hausmanniella discoidea, a type of the hausmanniellids (Bourland et al., 2011;Dunthorn et al., 2008).
The Paracolpoda genus was assigned to the Colpodidae and established by Lynn (1978), who noted that Colpoda species have a somatic groove while Paracolpoda species do not. Foissner (1985) subsumed C. steini into Paracolpoda. A close molecular relationship, highly similar merotelokinetal mode, somatic ciliary pattern and basic organization of the oral apparatus with P. steinii suggests Bromeliothrix metopoides could be temporarily assigned to Colpodidae (Foissner, 1993, 2010, Foissner et al., 2002, 2011. Additionally, the rather long, semicircular curved right oral polykinetid character, among other features of the Hausmanniellidae (Foissner, 1993) might be interpreted as convergence. Until now, the molecular data available in the NCBI database for Bromeliothrix and Paracolpoda are very sparse, with sequences from only one species for each. Broader sampling and molecular sequencing of Paracolpoda species is necessary to further clarify these relationships. Paracolpoda steinii pop5 separates from P. steinii pop1-4 in concatenated tree, the reason may be due to morphological differences, e.g., a larger body size compared with the other four populations.

| Bardeliella pulchra evolutionary position
In previous work, and in our study, B. pulchra was always the earliest divergent branch of Colpodida in the SSU-rRNA gene tree, which means it has been stable between Colpodida and Cyrtolophosidida (Bourland et al., 2011;Dunthorn et al., 2012;Rajter et al., 2020;Vďačný & Foissner, 2019). But, we established multiple gene trees and they all confirmed the instability of the position of B. pulchra in Colpodida.
Morphologically, three traits indicated that Bardeliella pulchra belongs to the Colpodida: (1) like most of colpodids, B. pulchra can divide in cysts; (2) B. pulchra's silverline pattern consists of highly ordered, comparatively large meshes extending between two ciliary rows each, and this silverline pattern is widely found in Colpodida species (Foissner, 1993;Foissner et al., 2011); (3) the oral ciliary fields of B. pulchra. Bardeliella and Colpoda are similar, i.e., the left oral ciliary field consisting of a shorter proximal portion with equidistantly spaced, monokinetidal rows .

CO N FLI C T O F I NTE R E S T
All authors declare that they have no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data presented in the study are deposited in the NCBI database (https://www.ncbi.nlm.nih.gov/) repository, accession numbers, lengths, and G&C contents are shown in Appendix Table A2.